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Abstract. The use of spectral methods to study birth-and-death processes was pioneered 
by S. Karlin and J. McGregor. Their expression for the transition probabilities was made ex- 
plicit by them in a few cases. Here we complete their analysis and indicate a few applications 
of their very powerful method. 
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1 Introduction 



1.1 Preliminaries 

In a paper overflowing with new as well as classical ideas, S. Karlin and J. McGregor [23] 
considered the birth-and-death process whose one step transition probability matrix P is given 
by 

/r po \ 

qi n pi o 

q2 r2 P2 

V '• '•■ '••/ 

where all indicated entries are nonnegative, po + r o < 1 an d Pi + T{ + q\ < 1, % > 1. These 
inequalities mean that there is a positive probability of jump to an (ignored) coffin state, or 
even distinct coffin states for each i such that the corresponding sum is not 1. 

The paper [23] acknowledges and builds on two different earlier approaches: one consists in 
exploiting the spectral theory of an operator in an appropriate Hilbert space, see [12, 21, 25] 
and the other one relies on a hands-on study of sample path properties, see [13, 18, 19]. 

They introduce the polynomials Qj{x) by the conditions Q_i(0) = 0, Qo(x) = 1 and using 
the notation 



Q(x) 



(Qo(x) 
Qi(x) 



they enforce the recursion relation 



FQ(x) = xQ(x) 
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to argue the existence of a unique measure tp(dx) supported in [—1, 1] such that 



Qi(x)Qj(x)ip(dx) = Sij. 



Then they obtain what is nowadays called the Karlin-McGregor representation formula 



7Ti 



x n Q i (x)Q j (x)ip(dx), 



where the constants irj are given by 



7T0 = 1, 



PoPi ■ •■Pj-i 
qiq 2 ■■■qj 



i>i- 



(1.1) 



Many general results can be obtained from the representation formula given above; as they 
remark, this is nothing but the spectral theorem. 

After establishing many such results the authors turn their attention to the study of the 
spectral properties in the case when one has essentially constant transition probabilities, more 
specifically 
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(1.2) 



Here the conditions on the parameters ro, po, q and p are 

ro > 0, po > 0, r +p <l, p,q>0, p + q=l. 

One usually puts qo = 1 — to — po and when qo > there is a loss of total mass of the initial 
probability distribution. 

If one denotes the spectral measure by ip(dx) and by m(z) its Stieltjes transform 



m(z) = J 



1 ^(dx) 



there is a classical method to obtain m(z), which we recall now. If n{z) denotes the Stieltjes 
transform of the measure going along with the "chopped matrix" obtained by deleting the first 
row and column of P then one proves that 

m(z) = -1/ (z -r + qpon(z)). 

The computation of n{z) itself can be made using the same trick; moreover since (in the case 
of an essentially constant transition matrix) the "twice chopped matrix" is the same as the 
"chopped matrix" we get that n(z) satisfies n(z) = — l/(z + qpn(z)) or n(z) 2 pq + zn(z) + 1 = 0, 
which gives, by solving this quadratic equation and observing that n(z) should vanish at infinity, 
n(z) = ( — z + \J z 2 — 4pqj / (2pq). 

Before proceeding, it is useful to note that the function n(z) is well defined as a single valued 
function in the complex plane from which we have removed the closed interval [— 2^/pq, 2^/pq]. 
As one approaches this cut from above, the function n(z) has a nontrivial imaginary part coming 
directly from the square root. This square root has positive values on the real axis to the right 
of the cut and negative values on the real axis to the left of the cut. This basic observation will 
play an important role at several points below. 
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Returning to the expression above one gets for m(z) 



m(z) 



1 




z + \J z 1 - 4pq) 



This expression may have zeros in the denominator giving rise to poles in m(z). Observe that 
none of these zeros can take place for values of z = x in the cut itself. This is, for instance, 
a consequence of the fact that along the cut the square root term gives a nontrivial imaginary 
contribution that cannot be cancelled by the other terms in the denominator which are real 
valued. 

A good way to study this function m(z) is to get rid of the square root in the denominator. 
In this fashion one gets the expression for m(z) already given in [23, Section 2ii], namely 



There is a classical way of recovering ip(dx) from m(z) which consists of considering m{z) as 
meromorphic function defined on an appropriate Riemann surface and inverting this Stieltjes 
transform. The measure consists of a nice density plus possibly some delta masses with certain 
weights. The function m(z) is well defined in the same region that makes n(z) single- valued. 
The continuous part of the measure results from looking at the values of the imaginary part 
of m(z) on one side of the cut. The discrete masses in ip(dx) come from the residues at the 
simple poles of m{z) but a more detailed analysis is given below. Clearly this inversion process 
requires a bit of careful computation. 

There are several references in the literature where the explicit computation of distributions 
associated to some special classes of Jacobi matrices is handled. One can mention [15] (see 
also [5]) for the case of periodic coefficients, [1, 2] for the case of parameters with constant 
limits, [9] for asymptotically periodic parameters and [6] in a different general context. 

For many other aspects of birth-and-death processes the reader can see [7, 11, 16, 20, 22, 24]. 
Many of these reference deal with continuous time. 

1.2 The story behind this paper 

In [23] one finds two examples done in some detail. Specifically they deal with the case when 
ro = 0, and the case when po = p, = q, which insures that the total mass of the initial 
probability distribution is preserved. For ro = and p = Po, which gives rise to a variant of the 
Wigner semicircle law, the authors also consider a two-periodic transition probability matrix. 

As preparation for a larger project we had planned to complete the analysis given in [23] by 
considering the case given in (1.2). Fortunately, during this process, Michael Anshelevich pointed 
us in the direction of some references in a different field, and after appropriate conversions it 
became clear that a full analyis of the problem posed in [23] was possible by using results in the 
literature. 

The most recent relevant paper is [26] in the context of "free probability" and the reader 
may also want to consult [4]. Our basic reference is an earlier paper [8] where Stieltjes inversion 
along the lines indicated above is avoided altogether. The method in [8] consists of using an old 
method of Christoffel which can be found in the very classical reference [28], and then getting 
a nice corollary which amounts to using the method of Christoffel in reverse. 

In [26] the authors observe that there is a small error in [8] amounting to having normalized 
the continuous and the discrete parts of ip(dx) in different fashions. The method used in [26] is 




(1.3) 
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the method of Stieltjes indicated above. In fairness one should say that in [26] there are several 
details left out but that the final result, essentially already given in [8], is correct. See also [3]. 

One should remark that neither of the references [8, 26] mention the results in [23]. A look 
at more recent papers indicates that the results in [23] are not widely known and that several 
of these cases are being rederived separately. For instance, [27] does not refer to any of the 
papers [8, 23, 26] although the author deals with several special cases of the problem considered 
in these references, such as the Wigner, Marchenko-Pastur and Kesten-McKay measures. 

Going back to the problem considered in [23] and by appropriately adapting the results 
in [8, 26] one can see that the general expression for the measure ip(dx) is given as follows. The 
continuous part of the measure is given by 



— \J Apq — x 2 
2tt v 

[p - po)x 2 - r (2p - p )x + prl + p\q 



dx, \x\ < 2^fpq. 



To this continuous part one needs to add possibly one or two mass points, and at this point the 
analysis in [8, 26] becomes relevant. 

In this paper we do not reproduce the arguments given in [8, 26] to obtain the orthogonality 
measure. Instead we display the results in a "reader friendly" way, we go through the steps 
showing why certain residues turn out to be positive, and more importantly we exploit the spec- 
tral method discussed here to derive some information about "quasi-stationary distributions". 
This goes beyond the analysis in [8, 23, 26] and combines this analysis with important work by 
E. van Doom and P. Schrijner, see for instance [29, 30]. Some of these results may be hard to 
obtain unless one uses spectral methods as in this paper and gives a rather nice way to complete 
the analysis started in [23] a long time ago. 

Since our expressions differ slightly from those in [8] and [26] we have checked carefully the 
internal consistency of our choices, including all normalizations. This required determining the 
eigenfunctions of the matrix (1.2). Eventually we found out that [7, p. 204] referring to earlier 
work of Geronimus [14] had already essentially written down these eigenfunctions as 



Qi(x s ' 



2{P0 - p) T 3 (x*) + 2 ^U j (x*) - 

Po Po Po V ? 



(1.4) 



where Tj and Uj are the Chebyshev polynomials of the first and second kind, and x* = x/ (2^/pq). 

When the measure is made up of the continuous part given above plus the possible con- 
tributions coming from the residues at the simple poles of m(z), the squared norms of these 
polynomials are as in (1.1): 

i PoP j ~ 1 ■ ^ , 

7T = 1, TTj = — , j > 1. 

qj 



2 Beyond Karlin and MacGregor 
2.1 A single pole 

The denominator of m(z) has a single root exactly when po = p, r$ ^ 0, in which case the root 
is 

xi = r +pq/r . (2.1) 

Here is a simple argument showing that the root has to be to the right of the interval 
[— 2^/pq, 2y/pq\. Indeed, we must have that r^+pq/rQ should be at least 2^/pq, and the difference 
of these two quantities can be rewritten as {-y/ro — \/pq/ro) 2 . 
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Figure 1. Case p = 0.2. 



Figure 2. Case p = 0.5. 



In the boundary case vq = ^Jpq the simple root is located at 2^/pq and it will turn out that 
we have no mass. 

For m(z) we get the expression 



m(z) 



1 z — 2ro — \J z 2 — Apq 
2r z- (r +pq/r ) 



We are interested in the residue of m(z) at (2.1) and this calls for a careful evaluation of the 
numerator at z = x\. 

Under the square root we have (ro + pq/ro) 2 — Apq = (ro — pq/r^) 2 . Since we are using the 
determination of the square root function which is positive on the real axis to the right of 2yjpq 
and we have shown that x\ > 2^/pq we have for the residue the value 



l ( . pq 

7, — I r o H 2r 

2r V r o 



ro 



pq 



This can be rewritten as (1/2) ( — (1 —pq/r^) — |1 — pq/r^l), which gives —(1 — pq/rfy^, where 
the subindex + denotes, as usual, the positive part. Since the denominator in the definition 
of m{z) is x — z instead of z — x the weight attached to this root is given by the negative of the 
residue, namely 



pq 



(2.2) 



o/ + 



Here and later it is important to note that the presence of zeroes in the denominator of m(z) 
does not guarantee the presence of a mass at such a zero, since the residue may vanish. But if 
this residue does not vanish then it is automatically positive as we have seen in the simple case 
analyzed above. In order to have a mass when po = p one should have 



p < 1/2 



and 



^0 G {Vm, l-p] 



(2.3) 



see Figs. 1, 2 and 3. 
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Figure 3. Case p — 0.85. 



2.2 Two roots 

In the case when the denominator has two different zeros at locations Xi, i = 1,2, an analysis 
of the residues along the lines of the simpler case above gives for the values of the weights the 
expressions 



ipI 



\Xi - r 



0-2 - 4q(p - po) \p\xi - r \ 
In particular, we see that to have a mass at x% we need 

r\ > 4q(p - p ) 
as well as 



(2.4) 



(2.5) 



(^-r ) 2 <^. 

p 



(2.6) 



We will describe, in terms of the values of the parameters ro, Po and p, the different possi- 
bilities according to the number of positive masses we may have. We will show below that we 
have just one mass exactly when 



p < — and 2p 



'-ro < Po < 2p + \l- r o, 



or 



p>- and p >2p-.-r , 
2 V a 



ro > yjpq 



(2.7) 
(2.8) 



see Figs. 1 and 3. Here, the mass point is always located on the right of the cut. 
In the same way, we have two masses exactly when 

1 [p 
p<- and po>2p+J-r , 

2 V Q 



(2.9) 



see in particular Fig. 1. In this case the two mass points are located at both sides of the cut in 
[-l,l]\[-2y/pq,2y/pq]. 
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Indeed, we write the explicit expression of each Xj, % = 1, 2: 

r (2p - po) ± Po\/ r l ~ ^liP ~ Po) , 
20 -Po) 

In order to have a mass at Xj, i = 1, 2, considering (2.6), one needs to have 

<2^E|p-p |- (2.10) 



r o ± V r l ~ 4 <?(p ~ Po) 



To analyze carefully inequality (2.10) one has to consider separately the cases po < p and po > p. 

Suppose that po < p. We will show that in this case it is possible to have a mass only at X2 
when 

p >2p- A fe (2.11) 
V Q 

holds true. Indeed, if we suppose that there is a positive mass at xi, from (2.10) we can write 



ro + \J 'rl - 4q(p - po) < 2^(p-p ), 
thus one should have 

2^(p-p ) > r , 

that contradicts (2.5). Indeed, considering the previous inequality and then (2.5) one would 
have 

^(p- m )*<(l-^)r S , 

which is a contradiction. 

Analogously, to have a positive mass at X2 we write for (2.10) 



or equivalently 

r - 2^=(P - Po) < \Jrl~ Aq{jp-p Q ). 

If (2.11) holds true then the expressions in both sides of the previous inequality are nonnegative 
so we can consider the squared inequality, which happens to be equivalent to (2.11) as well. We 
remark that when po < p (2.11) implies ro > y/pq. 

We conclude the analysis of the case po < p showing that if X2 has a positive mass then 
X2 > 2Jpq. Indeed, the inequality X2 — 2Jpq > is equivalent to 



r (2p-p ) - ^\/pq{p ~ Po) > Po\Jrl -4q(p-p ). 

Using (2.5) one verifies that the left part of the previous inequality is nonnegative thus taking 
the squares and simplifying conveniently one obtains the equivalent expression 

{roy/p ~ Vv( 2 P ~ Po)) 2 > °- 
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We point out that in the boundary case po = 2p — \Jp/qr$ the weight for %2 = 2^fpq in (2.4) is 
zero. 

For the case po > p, using similar arguments as before, one can see that to have a mass at x\ 
one needs to have 

p >2p+^r . (2.12) 
V 1 

Since 

Po + r <l, (2.13) 

condition (2.12) implies p < 1/2. Analogously, if (2.11) holds true, there is a positive mass 
at X2- Consequently, in this case one has two masses when (2.9) holds true, and one verifies 
that x\ < —2^/pq and X2 > 2^/pq. 

Notice that, in both cases described above, in order to have a mass one should have p ^ 
1/2, otherwise we will get a contradiction with condition (2.13) since the parameters must 
verify (2.11). See Figs. 1, 2 and 3. 

One can consider the case when our process is recurrent (see Sections 2 in [23, 29, 30]) for 
characterizations of some probabilistic properties for this kind of model) , which holds when 

r +p Q = l and p < \. (2.14) 

In this case z = 1 is a root of the denominator in m(z), and has a mass here only when the 
process is positive recurrent as noticed in [23, Section 2], which corresponds to p < 1/2 in (2.14). 
As a very special subcase, one considers tq = 0, po = 1, corresponding to a reflecting boundary 
condition at the left most state. Then one has roots of m(z) both at z = 1 and z = —1, having 
weights only when p < 1/2, as pointed out initially in [23, Section 2] and [22, pp. 108-109]. 



2.3 Some useful figures 

Here we have included some useful figures illustrating the different situations that may arise 
for the roots of the denominator of m(z) in (1.3), discussed above. The nature of the figures 
depends on the value of a fixed value of the parameter p being less, equal or larger than 1/2. We 
consider free nonnegative parameters po and ro satisfying (2.13). As noticed before, in the case 
of Fig. 1 for p < 1/2 one can see that there may be one or two mass points. These situations 
are delimited by the lines po = 2p — ^Jp/qro and po = 2p + \/p/qro . We point out that in order 
to have two masses it is necessary that po > p. 

For p = 1/2 Fig. 2 shows that we do not have any mass points. For p > 1/2, as Fig. 3 shows, 
there may be one mass point above the line po = 2p — \/p/qro and for po < p. 



3 Miscellaneous 

3.1 Quasi-stationary distribution 

A new notion enters at this point. A quasi-stationary distribution for a Markov chain with an 
absorbing state (the case of our chain when qo does not vanish) is defined as an initial distribution 
such that the (conditional) probability of the chain being at state j at time n, conditioned on 
the fact that no absorption has occurred by then, is independent of the time n for all states j. 

The issue of the existence and/or uniqueness of such distributions can be illuminated by the 
analysis given above. 
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Assume that go > so that we have a loss of total mass for any initial probability distribution. 
Assume further that p < 1/2 so that eventual transition out of the state space is guaranteed. 
We denote by p_i(n) the probability of having left the state space by time n. 

One defines a (normalized) probability distribution a« on the nonnegative integers as a quasi- 
stationary distribution for our walk if, by denoting with Pj(n) the probability that the walk be 
at state j at time n for the random walk with initial distribution Pj(0) = ay, one gets the relation 

Pi (n) 

-t.-, for j = 0,1,2,..., n = 0,1,2,.... 



l-p_i(n) 

Denote now by 77 the supremum of all points in the support of the orthogonality measure. 
This can be a point mass to the right of 2^/pq or the value I^Jpq itself. In both cases one can 
have 7] = 1. This value plays a very important role in terms of the existence and uniqueness of 
quasi-stationary distributions as follows from a more general result in [29, Theorem 4.2]. 

In the case when qo > and when absorption at —1 is certain, which holds for p < 1/2 
(see [29, Theorem 2.2]), one gets that if 77 = 1 then there is no quasi-stationary distribution. In 
our model this situation only occurs for p = 1/2. Indeed, when p < 1/2, po + ro < 1 and po ^ p, 
x = 1 is a root of the denominator of m(z) in (1.3) if 

Po 

1 + Po ~ r 

But in this case one never has a mass at x = 1, since the positive part of the corresponding 
residue vanishes in (2.4). When p < 1/2, po + r o < 1 an d p = po, x = 1 is a root of the 
denominator of m(z) when ro = po- In this case the value of the weight in (2.2) is equal to zero. 

In case that rj < 1, which occurs in our model when p < 1/2 and qo > 0, we have a one 
parameter family of such distributions parametrized by values of x in the interval rj < x < 1. 

The actual values of the distributions (which depend on x) are given by 

/ \ / -1 \Qj\ x ) 
a>jlx) = 7T 7 -(1 — x)— , 

qo 

where Qj(x) are the eigenfunctions given in (1.4). 
3.2 Ratio limits 

From the very beginning of the work of Karlin and McGregor, involving birth and death processes 
with discrete or continuous time parameter, one of the most important advantages of the spectral 
approach was the possibility to study the value of the limit as n grows of quantities such as 

P^/Pll- (3.1) 

As [23] points out these are difficult quantities to study and back in the 1930's Doeblin intro- 
duced tools to study the behaviour of the ratio of partial sums of these quantities. The work 
in [23] can thus be seen as getting a Tauberian theorem to go beyond Doeblin's results. In [23, 
Theorem 2] one finds a proof of the fact that for a recurrent (and aperiodic) walk the limit 
above is given by the ratio 71^/717. In our case, the walk is aperiodic or non-symmetric when 
ro 7^ and recurrent when (2.14) occurs. 

Assume that the parameters ro,po,p,ro 7^ 0, are such that we have at least one nonzero 
point mass. From Sections 2.1 and 2.2 we know that we have just one mass when either p = po 
and (2.3) holds or p 7^ po along with conditions (2.7) or (2.8) and one finds two masses under 
condition (2.9). Then one can use a result in [30, Theorem 3.1] to see that the limit exists for 
all values of the indices i, j, k, I and is equal to 

nj Qi{r])Qj{r]) 



717 Qk(v)Qi(v) ' 



(3.2) 
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where r\ is the supremum of all points in the support of the spectral measure, as defined above. 
One can also see that for any values of the parameters ro,Po,P, ro 7^ 0, the limits of the subse- 
quences resulting by taking n even and odd respectively in (3.1) do exist. For the symmetric 
case ro = 0, these limits where found in [23, Theorem 3] to be equal to (3.2) for i — j and k — I 
being simultaneously even or odd. 

3.3 The matrix valued case 

We conclude by calling attention to the last section in [23] where the authors consider an 
extension of their spectral method to the study of a walk on the integers. They show quite 
explicitly that the orthogonality measure is to be replaced by a matrix valued measure, which 
they write down. This point is noticed, for instance, in the Math. Reviews note on this paper, 
MR0100927. In [23] one does not find the corresponding matrix valued orthogonal polynomials 
or the extension to the matrix valued case of the Karlin-McGregor formula. This was eventually, 
and independently done much later in [10] and [17]. 
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